#include "./cell.hpp"
#include "yocto/string/vfs-utils.hpp"
#include "yocto/ios/ocstream.hpp"
#include "yocto/string/conv.hpp"
#include "yocto/code/utils.hpp"
#include "yocto/eta.hpp"
#include "yocto/duration.hpp"
#include "./manip.hpp"
#include "yocto/math/fcn/zfind.hpp"

#if 0
static void save_header( const Cell &cell, const string &filename )
{
    ios::ocstream fp( filename, false );
    fp("#t");
    fp(" Em");
    fp(" pH");
    fp(" P_CO2");
    cell.sol.save_header(fp);
    fp("\n");
}

static void save_values( const Cell &cell, const string &filename, double t )
{
    ios::ocstream fp( filename, true );
    fp("%.15g",t/60.0);
    fp(" %.15g", cell.Em);
    fp(" %.15g", cell.sol.pH());
    fp(" %.15g", cell.P_CO2(t) * 760);
    cell.sol.save_values(fp);
    fp("\n");
}


static inline
void progress_callback( const eta &ETA, const double t )
{
    static const char   ch[] = { '|', '/', '-', '\\' };
    static const size_t nc   = sizeof(ch)/sizeof(ch[0]);
    static size_t       ic   = 0;
    const duration      done( ETA.time_done );
    const duration      left( ETA.time_left );
    
    
    fprintf( stderr, "[%c]t=%9.1fs | %5.1f%% in %02dh%02dm%04.1fs : ETA   %02dh%02dm%04.1fs    \r",
            ch[ ++ic % nc ], t, ETA.ratio_done * 100.0,
            done.h,
            done.m,
            done.s,
            left.h,
            left.m,
            left.s
            );
    fflush( stderr );
}
#endif


//! pressure driver
class Pressure
{
public:
    double P0;
    double t0;
    double P1;
    double t1;
    double slope;
    Cell  *cell;
    double h1;   //!< target h
    vector<double> state;
    
    Pressure() : P0(0), t0(0), P1(0), t1(0), slope(0), cell(0), h1(0)
    {
    }
    
    ~Pressure() throw() {}
    
    void compute_slope()
    {
        slope = (P1-P0)/(t1-t0);
    }
    
    double getP( double t)
    {
        const double ans = P0 + (t-t0) * slope;
        //std::cerr << "P(" << t << ")=" << ans << std::endl;
        return ans;
    }
    
    void attach( Cell &a_cell )
    {
        cell = &a_cell;
        state.make( cell->nvar, 0.0);
    }
    
    double get_delta_h( double P )
    {
        assert(cell!=NULL);
        assert(state.size()>=cell->nvar);
        P1 = P;
        compute_slope();
        cell->save_state(state);
        
        cell->Step(t0, t1);
        const double cell_h1 = cell->sol["H+"];
        //std::cerr << "P1=" << P*760 << " => pH1=" << -log10(cell_h1) << " / " << -log10(h1) << std::endl;
        cell->load_state(state);
        return cell_h1 - h1;
    }
    
    
    static Pressure *Handle;
    static
    int P_CO2( lua_State *L )
    {
        assert(Handle!=NULL);
        const double t = lua_tonumber(L,1);
        lua_pushnumber(L,Handle->getP(t));
        return 1;
    }
    
private:
    YOCTO_DISABLE_COPY_AND_ASSIGN(Pressure);
};

Pressure * Pressure::Handle = 0;

// TODO: modify
static void save_header( const Cell &cell, const string &filename )
{
    ios::ocstream fp( filename, false );
    fp("#t");
    fp(" Em");
    fp(" pH");
    fp(" P_CO2");
    cell.sol.save_header(fp);
    fp("\n");
}

static void save_values( const Cell &cell, const string &filename, double t )
{
    ios::ocstream fp( filename, true );
    fp("%.15g",t/60.0);
    fp(" %.15g", cell.Em);
    fp(" %.15g", cell.sol.pH());
    fp(" %.15g", cell.P_CO2(t) * 760);
    cell.sol.save_values(fp);
    fp("\n");
}


int main( int argc, char *argv[] )
{
    const char *progname = _vfs:: get_base_name(argv[0]);
    try
    {
        
        ////////////////////////////////////////////////////////////////////////
        //
        // check arguments
        //
        ////////////////////////////////////////////////////////////////////////
        if( argc < 4 )
            throw exception("usage: %s config.lua input_file output_file", progname );
        const string config_file( argv[1] );
        const string input_file(argv[2]);
        const string output_file(argv[3]);
        
        ////////////////////////////////////////////////////////////////////////
        //
        // create data
        //
        ////////////////////////////////////////////////////////////////////////
        Cell  cell(config_file);
        Manip manip(input_file);
        
        
        std::cerr << "cell.sol: @pH=" << cell.sol.pH() << std::endl;
        std::cerr << cell.sol << std::endl;
        
        std::cerr << "cell.out: @pH=" << cell.outside(0).pH() << std::endl;
        std::cerr << cell.outside(0) << std::endl;
        
        ////////////////////////////////////////////////////////////////////////
        //
        // recompute initial pH from data
        //
        ////////////////////////////////////////////////////////////////////////
        {
            const double h0  = manip.h[1];
            initializer &ini = cell.InitIn;
            bool could_change = false;
            for( size_t i=1; i <= ini.size(); ++i )
            {
                constraint &cst = ini[i];
                if( 1 == cst.size() && cst.has("H+") )
                {
                    double &w = *cst.weight_address_of("H+");
                    std::cerr << "Changing Initializer Internal pH to " << manip.pH[1] << std::endl;
                    cst.set_value( h0 * w );
                    could_change = true;
                    break;
                }
            }
            if( !could_change )
                throw exception("couldn't change pH");
            
            cell.ReloadInside();
            std::cerr << "cell.sol: @pH=" << cell.sol.pH() << std::endl;
            std::cerr << cell.sol << std::endl;
        }
        
        ////////////////////////////////////////////////////////////////////////
        //
        // compute steady and/or set steady
        //
        ////////////////////////////////////////////////////////////////////////
        std::cerr << "ComputeSteadyEm" << std::endl;
        cell.ComputeSteadyEm();
        
        {
            lua_State *L = cell();
            lua_settop(L, 0);
            lua_getglobal(L, "Em");
            if( lua_isnumber(L, -1) )
            {
                const double Eini = lua_tonumber(L, -1);
                std::cerr << "Detect user's Em=" << Eini << std::endl;
                cell.EnforceSteadyEm( Eini );
                cell.RescaleEffectors();
            }
        }
        
        cell.TryRescaleAllRates();
        
        
        
        ////////////////////////////////////////////////////////////////////////
        //
        // Get initial pressure
        //
        ////////////////////////////////////////////////////////////////////////
        Pressure     pressure;
        Pressure::Handle = &pressure;
        pressure.attach( cell );
        
        pressure.P0  = cell.P_CO2(0.0);
        pressure.t0  = manip.t[1];
        
        std::cerr << "Starting P0=" << pressure.P0 * 760 << std::endl;
        
        ////////////////////////////////////////////////////////////////////////
        //
        // override P_CO2 in lua
        //
        ////////////////////////////////////////////////////////////////////////
        {
            lua_State *L = cell();
            lua_settop(L, 0);
            lua_pushcfunction(L, Pressure::P_CO2);
            lua_setglobal(L, "P_CO2");
        }
        
        ////////////////////////////////////////////////////////////////////////
        //
        // generic function to be zero
        //
        ////////////////////////////////////////////////////////////////////////
        numeric<double>::function delta_h( &pressure, & Pressure::get_delta_h );
        zfind<double>             zsolve( 1e-4 );
        
        ////////////////////////////////////////////////////////////////////////
        //
        // reconstruct
        //
        ////////////////////////////////////////////////////////////////////////
        const double P_grow   = 1.1;
        const double P_shrink = 0.9;
        
        
        save_header(cell, output_file);
        save_values(cell, output_file, manip.t[1]);
        
        for( size_t i=2; i<= manip.t.size(); ++i )
        {
            
            pressure.t1 = manip.t[i];
            pressure.t0 = manip.t[i-1];
            //std::cerr << "t: " << pressure.t0 << " -> " << pressure.t1 << std::endl;
            pressure.h1 = manip.h[i];
            double   P_lo = pressure.P0 * P_shrink;
            double   P_hi = pressure.P0 * P_grow;
            while( delta_h(P_lo) * delta_h(P_hi) >= 0.0 )
            {
                P_lo *= P_shrink;
                P_hi *= P_grow;
            }
            const double P1 = zsolve( delta_h, P_lo, P_hi );
            //const double dh = delta_h(P1);
            
            
            //------------------------------------------------------------------
            // forward cell state
            //------------------------------------------------------------------
            pressure.P1 = P1;
            pressure.compute_slope();
            cell.Step( pressure.t0, pressure.t1);
            
            //------------------------------------------------------------------
            // get ready for next step
            //------------------------------------------------------------------
            pressure.P0 = P1;
            std::cerr << "P(" << pressure.t1 << ")=" << P1 * 760 << " mmHg" << std::endl;
            save_values(cell, output_file, manip.t[i]);
            
        }
        
        
        
#if 0
        
        
        
        
        ////////////////////////////////////////////////////////////////////////
        //
        // prepare to run
        //
        ////////////////////////////////////////////////////////////////////////
        const size_t num_iter = max_of<size_t>(1,ceil(tmax/dt));
        double       t = 0;
        eta          ETA;
        const double iter_factor = 1.0 / num_iter;
        //----------------------------------------------------------------------
        // save at 0
        //----------------------------------------------------------------------
        save_header(cell, output );
        save_values(cell, output, t);
        
        //----------------------------------------------------------------------
        // run
        //----------------------------------------------------------------------
        std::cerr.flush();
        ETA.reset();
        double tmx_old = -1;
        for( unsigned iter=1; iter <= num_iter; ++iter )
        {
            const double t_old = (iter-1) * dt;
            t                  =  iter    * dt;
            cell.Step(t_old, t);
            save_values(cell, output, t);
            const double tmx_now = ETA.now();
            if( tmx_now - tmx_old >= 0.1 )
            {
                ETA(iter*iter_factor);
                progress_callback(ETA,t);
                tmx_old = tmx_now;
            }
        }
        ETA(1.0);
        progress_callback(ETA,t);
        fflush( stderr );
        std::cerr << std::endl;
#endif
        
    }
    catch( const exception &e )
    {
        std::cerr << e.what() << std::endl;
        std::cerr << e.when() << std::endl;
        return -1;
    }
    catch(...)
    {
        std::cerr << "Unhandled exception in " << progname << std::endl;
        return -1;
    }
    
    
    return 0;
}
